Heat conduction and Fourier's law 
by consecutive local mixing and thermalization 
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We present a first-principles study of heat conduction in a class of models which exhibit a new 
multi-step local thermalization mechanism which gives rise to Fourier's law. Local thermalization in 
our models occurs as the result of binary collisions among locally confined gas particles. We explore 
the conditions under which relaxation to local equilibrium, which involves no energy exchange, takes 
place on time scales shorter than that of the binary collisions which induce local thermalization. 
The role of this mechanism in multi-phase material systems such as aerogels is discussed. 
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The transport property of thermal conductivity has 
been introduced as a phenomenological law by Joseph 
Fourier in 1822 [lj. Its molecular origins were later dis- 
covered with Boltzmann's and Maxwell's works on the 
kinetic theory of gases [2j, |3|. In metallic solids, the ther- 
mal conduction is mainly determined by the transport 
of electrons and thus related to the electric conductivity 
according to the Wiedemann-Franz law, as explained by 
Drude in 1900 4]. In contrast, the thermal conductiv- 
ity of non-metallic solids results from the mutual scat- 
tering of sound waves due to the anharmonicities of the 
inter-atomic forces, as shown by Peierls in 1929 0, 
In all these systems, the microscopic mechanism of ther- 
mal conduction is the scattering of the energy carriers, 
which induces local thermalization in a single-stage pro- 
cess, and subsequent uniformization of temperature Q. 
Yet, this one-stage scenario does not exhaust the pos- 
sible microscopic mechanisms. Here, we prove that, in 
a further class of systems, local mixing can precede the 
start of energy transfer leading to local thermalization. 
On the basis of this mechanism, we provide a derivation 
of Fourier's law. 

This new mechanism of local thermalization allows us 
to understand how materials can become excellent ther- 
mal insulators, as is the case with aerogels 0. In these 
materials, a gas is trapped in a solid-state nanoporous 
matrix. The gas molecules collide with the walls of 
the nanopores more frequently than among themselves. 
Since the collisional transfer of energy is smaller with a 
stiff wall than with another free molecule, the changes 
of velocity orientation occur before a significant energy 
transfer. Accordingly, we envisage the mechanism where 
the wall collision frequency is shorter than the binary 
collision frequency. 

In this respect, thermal conduction should proceed 
over three well-separated time scales, r wa n -C Tbi nary <C 
T macro , which are : (i) the short time scale T wa n of the 
collisions on the solid walls, inducing the mixing of the 
particle positions and velocity orientations with negligi- 
ble transfer of particle kinetic energy and no mass trans- 



port; (ii) the intermediate time scale Tbinary of the bi- 
nary collisions between pairs of particles, achieving a lo- 
cal thermodynamic equilibrium at some locally defined 
temperature; and (iii) the long time scale T macro of the 
macroscopic relaxation of the Fourier modes. 

To demonstrate this mechanism, we consider in this 
Letter Hamiltonian dynamical systems in which hard-ball 
particles undergo elastic collisions either with immobile 
obstacles composing the solid matrix, or with each other. 
The obstacles form a lattice structure of cells, each con- 
fining a single mobile particle and preventing any mass 
transport. The motion of the particles is controlled by 
the geometry of the cells in such a way that the binary 
collisions only occur between nearest-neighboring mobile 
particles at a rate which can be switched off by shrinking 
the domain of motion within each cell. At the critical 
geometry where the binary collisions become impossible, 
the thermal conductivity vanishes with the binary colli- 
sion frequency, in a way which can be rigorously calcu- 
lated, as we show below. The remarkable result is that 
Fourier's law of heat transport can be established for this 
class of systems, based on the strong ergodic properties 
that they enjoy 

Under the assumption that binary collisions are seldom 
compared to wall collisions, the chaotic motion of indi- 
vidual particles within their cells induces a rapid decay 
of statistical correlations. As a consequence, the global 
multi-particle probability distribution of the system typ- 
ically reaches local equilibrium distributions at the ki- 
netic energy of each individual particle before energy 
exchanges proceed. This allows us to obtain a kinetic 
equation for the probability distribution of the local en- 
ergies which, under a Boltzmann-type assumption, 
reduces to a Boltzmann-Kac equation [ll| . In contrast 
to the Boltzmann equation which requires a further ap- 
proximation, our kinetic equation can be fully justified 
thanks to the mixing property of the dynamics and the 
separation of the time scales r wa n <C Tbinary 

To be specific, we consider a system such as depicted 
in Fig. [TJ It consists of a periodic two-dimensional array 
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FIG. 1: (Color online) Example of lattice billiard with rhom- 
bic tiling. The mobile disks are color-coded from blue to red 
with growing kinetic energies. See text for details. 



of discs of radii pf. In every cell between those discs, 
we place a single mobile disc of radius p m , large enough 
that the particles are confined to the regions between 
the discs. The sha pe of this region, which is a bounded 
dispersing billiard [12j , only depends on the sum of the 
two radii p = pi + p m , which is kept constant. 

Energy exchanges occur when two moving particles lo- 
cated in neighboring cells collide. Such binary collisions 
can take place provided the radii p m of the moving par- 
ticles is large enough. For a lattice with cells of sides I, 
binary collisions are thus possible above the critical value 



p c = y/p 2 - (7/2) 2 < p m < p 



(1) 



In this regime, the system is a thermal conductor. On 
the contrary, it becomes a thermal insulator for p m < p c . 
The thermal conductivity vanishes as p m ^ p c in an 
exactly calculable way. 

In order to establish this result, we start from first 
principles. The system is of Hamiltonian character and 
the time evolution of the global probability distribution 
of the positions and velocities of all the N particles is 
described by the pseudo-Liouville equation [13J : 
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where — v a • d Ta is the advection term expressing the 
changes of the probability distribution due to the free 
motion at the velocity v a of particle a inside its cell, the 
wall term K^ a - k ^ is the operator ruling the collisions of 
particle a on the fixed disc k in the corresponding cell, 
and _B( a ' & ) rules the binary collisions between the par- 
ticles a and b. Without these last terms, the pseudo- 
Liouville equation would describe the relaxation of the 
probability distribution to the product of local distribu- 
tions Il^li ^( £ a — TOl a/2) f° r & 11 the particles at their 



respective kinetic energies e Q . This relaxation is driven 
by the local mixing of the single-cell dispersing billiards, 
which takes place over the short time scale r wa n. The 
important point is that the probability distribution has 
not yet reached a local Maxwellian distribution. The lo- 
cal thermalization is the outcome of the binary collisions 
which take place on the intermediate time scale Tbinary 
This stage of the process of heat conduction is in general 
utterly difficult to master in a mathematically control- 
lable way. However, the fact that local mixing precedes 
local thermalization allows us to obtain a kinetic equa- 
tion for the probability distribution of the energies of the 
particles : 

r N 

P^ q \e u ...,e N ,t)= /J]dr dv x (3) 

^ a=l 

N 

pjv(ri,vi, . . . , r N , v N , t) Y[ S(e a - mv 2 a /2), 



where v a = ||v a ||. Integrating the pseudo-Liouville equa- 
tion {2} with the local equilibrium distributions, the ad- 
vective terms as well as the terms of the wall collisions 
vanish. We remain with the terms due to the binary col- 
lisions, which yield the following master equation [1 01 ] : 



n f 

d t P^(e 1 ,...,e N ,t)^^ £ I dr, 

a, 6=1 



(4) 
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(2tt)2 \C p>Pm (2)| 
/ d^dR / dv a dv b e ab ■ v ab 8 (e a - — u^J x 

6 ( 6b ~ y Ufc2 ) 6 ( v ~ y ^ ab ' Va)2 ~ ^ ab ' Vfo)2] ) ' 

is the transition rate of the binary collisions of energy 
transfer r\ between the neighboring particles a and b. In 
this expression, \C PtPin (2)\ is the volume of the billiard 
corresponding to two neighboring cells a and 6, which 
can be approximated by the square of the volume of a 
single cell \B P \ 2 , <j) is the angle of the unit vector e ab — 
(cos <j>, sin (j>) connecting a and 6, and R = {r a +T b )/2 the 
center of mass of the particles a and b. The transition rate 
W can be further expressed in terms of Jacobian ellip- 
tic functions 14j . The master equation ((4]) describes the 
time evolution as a continuous-time stochastic process of 
Poisson type with exponential probability distributions 
for the waiting times between the random events of en- 
ergy transfers due to the binary collisions. Accordingly, 
the master equation ^ shows that local thermalization 
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is reached over the intermediate time scale Tbmary — 1/^b 
of the binary collisions, which, in terms of the binary col- 
lision frequency, is given by 



J de Q de 6 A-q W {e 



k B T 2/9 m 
irm \B p \ 2 



d</>dR 



(6) 



where = l/(k-QT) is the inverse of the temperature T 
multiplied by Boltzmann's constant k B . 

Starting from the master equation ([J), we derive the 
heat equation 



d t T = d x ( K d x T), 



(7) 



for the local temperature T = T a defined in terms of 
the kinetic energy of particle a averaged over the prob- 
ability distribution ([3]), (e Q ) = k B T a . The derivation of 
Eq. |(7J) from stochastic models such as defined by the 
master equation @ requires the computation of corre- 
lation functions [la ]. In particular the heat current be- 
tween two neighboring cells involves an average with re- 
spect to the joint probability distribution of their energies 
and thus depends on their correlations. Assuming non- 
equilibrium boundary conditions, a systematic computa- 
tion of the stationary state can be performed through a 
cluster expansion. Single cell distributions, which solve 
the Boltzmann-Kac equation (llj , give the first order so- 
lutions of this expansion, while pair correlations appear 
at the second order, and so on. For the stochastic system 
described by Eq. (j4|), it turns out that pair correlations 
are of second degree in the local temperature gradients 
and can therefore be neglected in the computation of the 
transport coefficient, which involves only terms linear in 
this quantity. We can thus limit the cluster expansion 
to first order and find the expression of the coefficient of 
thermal conductivity : 
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(8) 



We verified this result by direct numerical simulations of 
Eq. Q . Alternatively, n can be defined through a Green- 
Kubo formula. The result (jHJ) shows that k is given by 
the S(t) part of the energy current-current correlation dis- 
cussed in Ref. [15| . We notice that this result holds here 



even though the heat current does not have the gradient 
form of [la], i. e. it cannot be written as the difference 
of two local functions. The details will be presented else- 
where (l6| . 

The comparison between Eqs. © and ([5]) shows the 
equality of the thermal conductivity and collision fre- 
quency of the stochastic process described by Eq. ([4]). 
By extension, the same holds of the corresponding prop- 
erties of the billiard system in the limit p m — > p c : 
lim Pm ^ Pc k/(1 i/b) = 1. Furthermore, the thermal con- 
ductivity can be explicitly calculated in this limit and 



shown to vanish as the third power of the difference 
Pm — Pc- Indeed, the geometric integral appearing in 
both Eqs. ((6]) and (|8|) vanishes as 



J d</>dR = c 3 (p n 



Pc) 3 + Ci(p m - p c ) 4 



(9) 



where the coefficients in this expansion can be computed 
according to the specific geometry of collision events [3] • 
Numerical computations of the quantities computed in 
Eqs. © and © are shown in Fig. [21 asserting the validity 
of the stochastic description on the one hand and the 
computation of its properties on the other. 




FIG. 2: (Color online) Thermal conductivity k in reduced 
units, computed by molecular dynamics from the heat ex- 
change of a single particle with stochastically thermalized 
neighbors, and binary collision frequency v^, as functions of 
p m — p c , compared to our theoretical calculations (J6]l and l|8|l 
(solid line). The inset shows the ratios between v\, and k and 
the RHS of ([6]) and p]). The computations were performed 
with a rhombic cell in contact with stochastically thermal- 
ized cells on each sides, in which case cz = 128p c /(3Z 2 ) and 
C4 = 256pc/(3^ 4 ) in Eq. (J9j . The sum of radii is here fixed to 
the value p — pt + p m = 9/25. 



Finally, the local mixing property can be character- 
ized by the Lyapunov exponents measuring the stretch- 
ing and contraction of phase-space volumes. Remarkably, 
all these exponents can be analytically calculated in the 
same critical limit where the binary collision frequency 
vanishes. In this limit, the particle velocities share a 
Maxwellian probability distribution at the inverse tem- 
perature (3 while the particles become independent. Ac- 
cordingly, a system containing TV particles presents the 

positive Lyapunov exponents 



In- 



N 



m/3 £ — 1/2 



(10) 



where A + is the positive Lyapunov exponent of a single 
isolated particle moving at unit speed. 2N exponents are 
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zero in the said limit and the N remaining exponents take 
the negative values — A,-. The sum of all the exponents 
vanishes since the system is Hamiltonian and satisfies 
Liouville's theorem. The system is thus chaotic with an 
extensive Kolmogorov-Sinai entropy per unit time given 
by h KS =N\ + y/Tr/(2mt3). 

Figure [3] shows the agreement between the numerically 
computed Lyapunov exponents and Eq. (|10|) . in the limit 
where the system becomes insulating. 
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FIG. 3: (Color online) Lyapunov exponents Xi versus their 
index i. The numbers in the legend correspond to the val- 
ues of p m . The first one, p m = 0.04 < p c , yields the pos- 
itive exponent A+ associated to isolated cells, with all the 
particles at the same speed. The solid red line is the first 
half of the spectrum as predicted by Eq. QD) . The inset 
shows Xi versus p m — p c for the first half of the positive expo- 
nents, Ai, . . . , Ajv and compares them to the asymptotic value 
(|10l) (dashed lines). The computations were performed with 
a rhombic channel of size TV = 1 cells, p = 9/25. 

To conclude, we have presented in this letter a mecha- 
nism of heat conduction which proceeds in two stages 
and involves no mass transport. The first stage is a 
local mixing phase which takes place in the individual 
cells and involves no energy transfer between them. This 
stage thus induces the relaxation of phase space distri- 
butions to local equilibrium without thermalization. The 
second stage is the energy transfer which leads to local 
thermalization and can be efficiently described in prob- 
abilistic terms by means of our master equation. The 
uniformization of the temperature over the whole sys- 
tem can then proceed as described by the heat equation. 
Whereas the first-principles derivation of Fourier's law is 
problematic based on the standard one-stage local ther- 
malization mechanism, it is remarkable that it can be 
carried out in the presence of this additional local equilib- 
rium mechanism. We notice that both mechanisms cor- 
respond to physically distinct limits for the heat conduc- 
tivity. In the standard one-stage mechanism, the ballistic 
transport of the energy carriers is hampered by scatter- 
ing so that the theory should there demonstrate a finite 
conductivity with respect to an unperturbed situation 



where the conductivity is unbounded. In contrast, in the 
two-stage mechanism we here describe, the thermal con- 
ductivity vanishes in the reference system, allowing the 
first-principles derivation of Fourier's law. This two-stage 
mechanism is encountered in aerogels where the thermal 
conductivity takes the smallest known values. Our anal- 
ysis therefore suggests that fundamental insight into the 
phenomenon of heat conductivity could be obtained by 
the experimental study of aerogels. 
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